# coding=utf-8
#version  liantongyu-youhua_0315_2210.py

import numpy as np
import cv2
from matplotlib import pyplot as plt
# from graythresh_own import graythresh_own
from skimage import io,data,color,measure,filters,draw,morphology
from skimage.measure import label, regionprops, regionprops_table
import math
from scipy.signal import convolve2d
import time
import sys
from time import gmtime, strftime
# from skimage.color import rgb2gra

def liantongyu2(kmeansFilename,choseSec=1):

    # img0=plt.imread('input.jpg') #读取原图像
    # img0Gray = Image.open('input.jpg').convert('LA')

    # img22=io.imread(kmeansFilename)
    # img11=io.imread(originimg)
    # tempImg=io.imread(kmeansFilename,True)
    # Gray = R*0.299 + G*0.587 + B*0.114
    # tempImgHH=tempImg[:,:,0]*38+tempImg[:,:,1]*75+tempImg[:,:,2]*15
    # BWimg=np.right_shift(tempImgHH,7)
    # BWimg=tempImg[:,:,0]/255*0.299+tempImg[:,:,1]/255*587+tempImg[:,:,2]/255*0.114

    # tempImg=io.imread(kmeansFilename,True)
    # plt.imshow(tempImg)
    # plt.show()
    # plt.imshow(BWimg)
    # plt.show()
    BWimg=io.imread(kmeansFilename,True)
    # BWimg = rgb2gray(BWimg)

    # # 基于Otsu的阈值分割方法
    thresh = filters.threshold_otsu(BWimg)   #返回一个阈值
    BWimg =(BWimg <= thresh)*1.0   #根据阈值进行分割

    # # %先闭运算 再开运算
    # se=strel('disk',2);
    # BWimg = imclose(BWimg,se);
    # BWimg = imopen(BWimg,se);
    # BWimg=morphology.binary_closing(BWimg,morphology.disk(2))  #用边长为5的圆形滤波器进行膨胀滤波
    BWimg=morphology.binary_opening(BWimg,morphology.disk(4))  #用边长为5的圆形滤波器进行膨胀滤波
    # cross = np.array([[0,1,0],[1,1,1],[0,1,0]])
    # BWimg = morphology.binary_erosion(BWimg, morphology.square(2))

    BWimg=np.array(BWimg).astype(int) #将BWimg float64转换成int32

    label_img = label(BWimg,connectivity=1)
    regions = regionprops(label_img)
    rowCnts,colCnts=np.shape(label_img)
    label_total = max(np.reshape(label_img,rowCnts*colCnts)).astype(int)


    # for props in regions:
    #     # y0, x0 = props.centroid
    #     # orientation = props.orientation
    #     # x1 = x0 + math.cos(orientation) * 0.5 * props.minor_axis_length
    #     # y1 = y0 - math.sin(orientation) * 0.5 * props.minor_axis_length
    #     # x2 = x0 - math.sin(orientation) * 0.5 * props.major_axis_length
    #     # y2 = y0 - math.cos(orientation) * 0.5 * props.major_axis_length

    #     # ax.plot((x0, x1), (y0, y1), '-r', linewidth=0.3)
    #     # ax.plot((x0, x2), (y0, y2), '-r', linewidth=0.3)
    #     # props.label
    #     # ax.text(x0, y0,props.label,fontsize=6,color='red')
    #     # ax.plot(x0, y0, '.r', markersize=1.5)

    #     minr, minc, maxr, maxc = props.bbox
    #     # rec= props.bbox
    #     # print(rec)
    #     # ax.plot(minr, minc, '-r', linewidth=0.5)

    #     # # 显示被注释
    #     # plt.gca().add_patch(plt.Rectangle(xy=(minc-1, minr-1),
    #     # width=maxc-minc+1,
    #     # height=maxr-minr+1,
    #     # edgecolor='red',
    #     # fill=False, linewidth=0.5))
    #     # plt.text(x0, y0,props.label,fontsize=6,color='red')

    # choseSec=92
    choseImg=(label_img==choseSec)*1   # choseImg原来为bool类型，当*1转换后变为int32类型
    # choseImg=(label_img==choseSec)   # choseImg为bool类型
    # print(np.dtype(choseImg[1,2]))   # 打印choseImg元素类型

    minr, minc, maxr, maxc = regions[choseSec-1].bbox  # regions为所有区块的属性列表List类型
    # print(regions[choseSec-1].bbox)
    choseSecImg=choseImg[minr:maxr,minc:maxc]

    rowColOrigin=np.array([minr, minc])

    # # 显示被注释
    # plt.show()

    coreMatrix=np.array([[0,1,0],[1,0,1],[0,1,0]])
    ww = convolve2d(choseSecImg,coreMatrix,mode='same',boundary='fill',fillvalue=0)
    choseSecImgBorder=((ww<4) & (ww>0))*1

    # **********功能：获得边界的像素点索引值indexM
    indexM = np.array(np.where(choseSecImgBorder == 1))
    indexM = np.transpose(indexM)

    return rowColOrigin,label_total,choseSecImg,choseSecImgBorder,indexM

def getImgBorderErZhiHua(filename):
    print( "--- liantongyu , getImgBorderErZhiHua  %s--start1-" % ( filename ) )
    BWimg=io.imread(filename,True)
    #BWimg=cv2.imread(filename)
    #BWimg=cv2.imread(filename,cv2.IMREAD_GRAYSCALE)
    print( "--- liantongyu , getImgBorderErZhiHua  %s--start1-" % ( filename ) )
    # # 基于Otsu的阈值分割方法
    thresh = filters.threshold_otsu(BWimg)   #返回一个阈值
    BWimg =(BWimg <= thresh)*1.0   #根据阈值进行分割

    # # %先闭运算 再开运算
    # se=strel('disk',2);
    # BWimg = imclose(BWimg,se);
    # BWimg = imopen(BWimg,se);
    # BWimg=morphology.binary_closing(BWimg,morphology.disk(2))  #用边长为5的圆形滤波器进行膨胀滤波
    BWimg=morphology.binary_opening(BWimg,morphology.disk(4))  #用边长为5的圆形滤波器进行膨胀滤波
    # cross = np.array([[0,1,0],[1,1,1],[0,1,0]])
    # BWimg = morphology.binary_erosion(BWimg, morphology.square(2))

    BWimg=np.array(BWimg).astype(int) #将BWimg float64转换成int32

    coreMatrix=np.array([[0,1,0],[1,0,1],[0,1,0]])
    ww = convolve2d(BWimg,coreMatrix,mode='same',boundary='fill',fillvalue=0)
    BWimgBorder=((ww<4) & (ww>0))*1

    # # **********功能：获得边界的像素点索引值indexM
    # indexM = np.array(np.where(BWimgBorder == 1))
    # indexM = np.transpose(indexM)
    return BWimgBorder

def judgeInSection_Bool_inSec(imgsectlogic,x0,y0,x1,y1):
    # minusCnt=4
    inSec=True
    imgsectlogicCopy=np.copy(imgsectlogic)
    linePosIndex=draw.line(x0,y0,x1,y1)
    imgsectlogicCopy[linePosIndex]=1
    imgsectdiff=imgsectlogic-imgsectlogicCopy
    minusCnt=np.sum(imgsectdiff<=-1)
    if(minusCnt > 3):
        inSec=False
    return inSec

def findMaxDisInSec_Vec_XYPointMaxDis(xyDisPointM, choseSecImg, calVerti=False):

    # *********功能：升序排列所有距离值xyPointM，并找出在区域内部choseSecImg的线段
    xyDisPointMSort=xyDisPointM[np.argsort(xyDisPointM[:,4]),:]
    [rowsCnt,colsCnt]=xyDisPointMSort.shape
    maxPointIndex=0
    for i in range(1,rowsCnt):
        [y0,x0,y1,x1]=xyDisPointMSort[-i,0:4].astype(int)
        inSection=judgeInSection_Bool_inSec(choseSecImg,x0,y0,x1,y1)
        if(inSection):
            maxPointIndex=-i
            break
    [y0Max,x0Max,y1Max,x1Max]=(xyDisPointMSort[maxPointIndex,0:4]).astype(int)
    if (calVerti == True):
        disMax = (xyDisPointMSort[maxPointIndex, 4]).astype(float)
        sigmaAngle = (xyDisPointMSort[maxPointIndex, 5]).astype(float)
        pointAndDisInSecM = [x0Max, y0Max, x1Max, y1Max, disMax, sigmaAngle]
    else:
        disMax = (xyDisPointMSort[maxPointIndex, 4]).astype(float)
        pointAndDisInSecM = [x0Max, y0Max, x1Max, y1Max, disMax]
    return pointAndDisInSecM

def calDis_Mat_XYPointDis(indexM):

    # *******功能：计算所有点（数量n=rows）之间的欧式距离，totals=n(n-1)/2，存储于disXYM=[x0,y0,x1,y1,distemp]
    rows,cols=np.shape(indexM)
    totals=int(rows*(rows-1)/2)
    xyPointDisM=np.zeros((totals,5),float)
    hh=0
    for i in range(0,rows):
        [y0,x0]=(indexM[i,0:2]).astype(float)
        for j in range(i+1,rows):
            [y1,x1]=(indexM[j,0:2]).astype(float)
            distemp=((x0-x1)**2+(y0-y1)**2)**0.5
            xyPointDisM[hh,:]=[x0,y0,x1,y1,distemp]
            hh=hh+1
    return xyPointDisM

def calVerticalDis_Mat_XYPointDisAngle(indexM,maxPointAndDisInSecM):

    rows,cols=np.shape(indexM)
    # [y0Max,x0Max,y1Max,x1Max]=(maxPointAndDisInSecM[0:4])
    [x0Max,y0Max,x1Max,y1Max]=(maxPointAndDisInSecM[0:4])
    # OA为内部最长线段2点的向量，用于寻找垂直2点，如OB
    OA=np.array([y0Max-y1Max,x0Max-x1Max],float)
    # OA=np.array([x0Max-x1Max,y0Max-y1Max],float)

    # **********功能：计算所有垂直线段最长距离
    # dissigmaXY=[x0,y0,x1,y1,distemp,sigmatemp]
    dissigmaXY=np.zeros((rows-1,6),float)
    # xyVertical为找出的任意2点垂直线段的距离
    xyVertical=np.zeros((rows,6),float)

    for i in range(0, rows):
        [y0, x0] = (indexM[i, 0:2]).astype(float)
        gg = 0
        for j in range(0, rows):
            if (i != j):
                [y1, x1] = (indexM[j, 0:2]).astype(float)
                OB = np.array([x0 - x1, y0 - y1], float)
                Lx = (OA.dot(OA))**0.5
                Ly = (OB.dot(OB))**0.5
                #相当于勾股定理，求得斜线的长度
                cos_angle = OA.dot(OB) / (Lx * Ly)
                #求得cos_sita的值再反过来计算，绝对长度乘以cos角度为矢量长度，初中知识。。
                sigmatemp = math.fabs(cos_angle)  #用cos_angle直接替代sigmatemp会快一点
                # if(cos_angle>1):
                #     cos_angle=1
                # if(cos_angle<-1):
                #     cos_angle=-1
                # angle=math.acos(cos_angle)*180.0/np.pi
                # sigmatemp=math.fabs(angle-90)
                distemp = ((x0 - x1)**2 + (y0 - y1)**2)**0.5
                dissigmaXY[gg, :] = [x0, y0, x1, y1, distemp, sigmatemp]
                gg = gg + 1
        sigmaM = dissigmaXY[:, 5]
        ddVerIndex = np.argmin(sigmaM)
        xyVertical[i, :] = dissigmaXY[ddVerIndex, :]
    return xyVertical

# def drawLinesInImage(imageMat,pointVec,figurename='figure',colorVar='r'):
#     (r0, c0, r1, c1) = pointVec[0:4]
#     xm = [c0, c1]
#     ym = [r0, r1]
#     floatimg = np.array(np.copy(imageMat)).astype(float)
#     floatimg = floatimg * 1
#     plt.title(figurename)
#     plt.imshow(floatimg)
#     plt.plot(xm, ym, linestyle='-', linewidth=1, color=colorVar)
#     # plt.show()
# return

# def drawLinesPixInImage(imageMat,pointVec,figurename='figure',colorVar=0.7):
#     (r0, c0, r1, c1) = pointVec[0:4]
#     linePosIndex = draw.line(r0, c0, r1, c1)
#     # imageMatCopy = np.copy(imageMat)
#     floatimg = np.array(imageMat).astype(float)
#     floatimg = floatimg * 1
#     floatimg[linePosIndex] = colorVar
#     plt.title(figurename)
#     plt.imshow(floatimg)
#     # plt.show()
# return


def genRadiusImg( kmeansFilename, maxLabelNum, savefigPath ):
    print("--- liantongyu, start %s  %d   %s ---" % (kmeansFilename, maxLabelNum, savefigPath) )
    pointVecList=[]
    pointVecList.clear()
    # pointVecVerticalList=[]
    np.set_printoptions( threshold=np.inf )
    #imageMat=io.imread( kmeansFilename )
    BWimgBorder=getImgBorderErZhiHua( kmeansFilename )
    rowColOrigin,label_total,choseSecImg,choseSecImgBorder,indexM=liantongyu2( kmeansFilename,1 )
    for i in range(1,maxLabelNum):
        print( "--- liantongyu , %s  %s--current label %d and totol %d-" % ( sys._getframe().f_lineno, strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime()) , i, label_total) )
        rowColOrigin,label_total,choseSecImg,choseSecImgBorder,indexM=liantongyu2(kmeansFilename,i)

        # pointVec = np.array([10, 20, 100, 30])
        tt1=time.time()
        xyPointDisM = calDis_Mat_XYPointDis(indexM)
        inter01=time.time()-tt1

        tt1=time.time()
        pointVec = findMaxDisInSec_Vec_XYPointMaxDis(xyPointDisM, choseSecImg)
        inter02=time.time()-tt1

        tt1=time.time()
        xyVertical = calVerticalDis_Mat_XYPointDisAngle(indexM, pointVec)
        inter03=time.time()-tt1

        tt1=time.time()
        pointVec2 = findMaxDisInSec_Vec_XYPointMaxDis(xyVertical, choseSecImg,calVerti=True)
        inter04=time.time()-tt1

        inter11=[inter01,inter02,inter03,inter04]
        print(inter11)

        rowColOrigin=np.array(rowColOrigin).astype(int)
        pointVec=np.array(pointVec)
        pointVec2=np.array(pointVec2)

        r0=pointVec[0]+rowColOrigin[0]
        c0=pointVec[1]+rowColOrigin[1]
        r1=pointVec[2]+rowColOrigin[0]
        c1=pointVec[3]+rowColOrigin[1]
        dis0=pointVec[4]
        rv0=pointVec2[0]+rowColOrigin[0]
        cv0=pointVec2[1]+rowColOrigin[1]
        rv1=pointVec2[2]+rowColOrigin[0]
        cv1=pointVec2[3]+rowColOrigin[1]
        disv0=pointVec2[4]
        pointMaxDisAndVerDis=np.array([r0,c0,r1,c1,dis0,rv0,cv0,rv1,cv1,disv0])  #第一组参数为最长线坐标和距离，第二组为垂直最长线坐标和距离

        (r0, c0, r1, c1) = pointMaxDisAndVerDis[0:4]
        xm = [c0, c1]
        ym = [r0, r1]
        (r0, c0, r1, c1) = pointMaxDisAndVerDis[5:9]
        xv = [c0, c1]
        yv = [r0, r1]
        plt.plot(xm, ym, linestyle='-', linewidth=1, color='r')
        plt.plot(xv, yv, linestyle='-', linewidth=1, color='g')

        pointVecList.append(pointMaxDisAndVerDis)

    print(pointVecList)
    # np.savetxt('test.txt',pointVecList)
    plt.imshow(BWimgBorder)
    plt.savefig( savefigPath )
    #io.imsave('BWimgBorder.jpg',BWimgBorder)
    #plt.show()
    return pointVecList


#genRadiusImg('result2.jpg',20,'./test1.jpg')